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Abstract 

The major goal of this paper is to find accurate solutions for linear fractional differential equations 
of order 1 < a < 2 . Hence, it is necessary to carry out this goal by preparing a new method called 
Fractional Finite Difference Method (FFDM). However, this method depends on several important 
topics and definitions such as Caputo’s definition as a definition of fractional derivative, Finite 
Difference Formulas in three types (Forward, Central and Backward) for approximating the second 
and third derivatives and Composite Trapezoidal Rule for approximating the integral term in the 
Caputo’s definition. In this paper, the numerical solutions of linear fractional differential equations 
using FFDM will be discussed and illustrated. The purposed problem is to construct a method 
to find accurate approximate solutions for linear fractional differential equations. The efficiency of 
FFDM will be illustrated by solving some problems of linear fractional differential equations of order 
1 < a < 2. ©2016 All rights reserved. 

Keywords: Finite difference formulas, composite trapezoidal rule, numerical solutions, linear 
fractional differential equation. 
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1. Introduction 

Fractional calculus is an old mathematical concept dating back to 17th century and involves 
integration and differentiation of arbitrary order. In a letter on 30th September 1695, L’Hospital 
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wrote to Leibnitz asking him about the differentiation of order 1/2. Leibnitz response was ’’apparent 
paradox from which one day useful consequences will be drawn”. In the following centuries, fractional 
calculus developed significantly within pure mathematics |21j . However the applications of fractional 
calculus just emerged in last few decades in various areas of physics and engineering, namely in 
signal processing, control engineering EM, electromagnetism [I0], bioscience ra, fluid mechanics 
[16], electrochemistry EES. diffusion processes [28], dynamic of viscoelastic materials H , continuum 
and statistical mechanics [16] and propagation of spherical flames [18] . Motivated by increasing 
number of applications of fractional differential equations, considerable attention has been given to 
provide efficient methods for exact and numerical solutions of fractional differential equations. 

In general, most of the fractional differential equations do not have exact solutions. Particularly, 
there is no known method for solving fractional differential equations exactly. Therefore several 
methods for approximate solutions to classical differential equations are extended to solve differential 
equations of fractional order numerically. Some of these methods include Adomian Decomposition 
Method [7] , Perturbation Method p[] , Homotopy Analysis Method m , Variational Iteration Method 
ra, Extrapolation Method [8], and Generalized Differential Transform Method [19], 

Fractional Finite Difference Method is useful to solve the linear fractional differential equations 
of order 1 < a < 2 numerically. This method is a technique by which fractional derivatives of y(t ) 
are approximated by differences in the values of the y{t) between a value of the independent variable 
t and small increment ( t + h). To solve fractional differential equations numerically, we can replace 
the derivatives in the linear fractional differential equation of order 1 < a < 2 with finite difference 
approximations. 

Further, the FFDM can be used to fold accurate approximate solutions to the general fractional 
differential equation of the form 

D°y(t) + a m y {m \t ) + a m _ 1 y < ' m - 1) (t) + ... + a 0 y(t) + N(y(t),y'(t)) = f(t ), ^ ^ 

t > 0, m — 1 < a < m, ' ' 

Subject to 0) = y % , i — 0,1,..., m — 1. 

Where D 2 is the derivative of y of order a in the sense of Caputo fractional differential operator, y(t) 
is an unknown function of the independent variable t and N is a nonlinear differential operator. 

FFDM’s efficiency is demonstrated by comparing the obtained approximate solutions for some 
linear fractional differential equations of order 1 < a < 2 by that method with the exact solutions 
and with other approximate solutions obtained by other solvers and methods. For the linear frac¬ 
tional differential equations, the mentioned method is utilized to fold approximate solutions for these 
equations. Moreover, Mathematica (Version. 9.0) software program is used to solve the system of 
linear equations that is obtained from some operations. 

2. Numerical Approximations and Fractional Calculus 

The differentiation and integration is one of the most important topics in mathematics. This topic 
deals with the Integer Calculus (A calculus of the integer order derivatives and integrals). Integer 
Calculus is tightly related to the well-known fundamental theorem of classical calculus [231 EH- 
Thus, the obtained results of integer order derivatives and integrals can often be carried over to 
the fractional case [23] . Fractional Calculus is a branch of mathematical analysis that studies the 
possibility of taking real or complex powers as orders. 

The fractional derivative of order a > 0 has several definitions. Riemann-Liouville and Caputo’s 
definitions are the most commonly used for the derivative of this order. For the fractional derivative, 
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the Caputo’s definition is used, which is a modification of the Riemann-Liouville definition; because 
it has an advantage of dealing properly with the initial value problem since the initial condition is 
given in terms of the field variables and their integer order. This case is widely used in physical ap¬ 
plications [6j. However, some basic concepts, definitions and properties of numerical approximations 
and fractional calculus are presented to be used in this section. 


2.1. Numerical Differentiation and Integration 

Usually, Taylor Series is used to obtain numerical differentiation methods. Methods using Taylor 
Series are: Backward Difference Method, Forward Difference Method and Central Difference Method 
to evaluate the derivative Emu. Finite difference formulas are the most common formulas which 
are used to solve the ordinary and partial differential equations numerically H3- The derivatives 
in these equations can be replaced with suitable finite difference approximations on a discretized 
domain. The accuracy of the solution depends on the number of mesh points such that if the number 
of mesh points is increased, then the solution will be more accurate [12]. However, Taylor Series can 
be used to derive some of finite difference approximations. For / (x 0 ), the well known formula is: 


f " / ^ f ( x o + 2/i) — 2/ (x 0 + h) + f(x 0 ) 

f ( x o) = -- + 0{h). 

The central finite difference formula for the 3 rd derivative is known as 

/" , ^ —/ (a?o - 2 h) + 2/ (x 0 -h) -2f (x 0 + h) + f (x 0 + 2 h) , 

/ —-^73- r )■ 


( 2 . 1 ) 


( 2 . 2 ) 


The forward finite difference formula for the 3 rd derivative similarly, which can be expressed as 


j'" ^ _ ~/ (go) + 3/ (xp + h) - 3/ (x 0 + 2 h) + f (x 0 + 3 h) +Q ^ 


(2.3) 


In addition, if we set —h instead of h in Eq. 
the 3 rd derivative will be obtained as 


(2.3), then the backward finite difference formula for 


t»> f \ f (®o) - 3/ (x 0 - h) + 3/ (x 0 - 2h) - f (x 0 - 3h) 

/ (®o) =--+ 0(h) 


(2.4) 


On the other hand, there are several numerical methods for estimating definite integrals; one of these 
methods is called Composite Trapezoidal Rule [5], 


Theorem 2.1. Suppose that the interval [ a,b\ is subdivided into n subintervals [xk,Xk+i\ of width 
h = Of by using the equally spaced nodes Xk = a + kh , for k — 0, 1 , ..., n. The Composite 
Trapezoidal Rule for n subintervals can be expressed as 

h n 

T (/> h ) = 2 U + f fa*)] • 

z k= 1 

This is an approximation to the integral of f(x) over [a, b] since 

/ (x) Ax k T (/, h). 



(2.5) 
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2.2. Fractional Derivatives 

We start this subsection by illustrating some notations that will be used frequently. The notation 
D"/ ( t ) will be introduced as a notation which is denoted to the fractional derivative of a function 
/ ( t ) along t- axis of an arbitrary order a > 0, where the subscript a is denoted to the lower limit 
of the integration. However, to get a simplified notation, we may drop the subscript a. Further, 
the fractional derivative can be defined using the definition of the fractional integral. The following 
definition illustrates that Ea- 

Definition 2.2. Let a G M + . For a positive integer m such that m — 1 < a < m, then the 
Riemann-Liouville Fractional Derivative of a function f{t) of order a is defined by 

D ‘ J {t) = TJ^ha) (!) L {t - Xr ^' f {X)AX ■ 

In addition, one observation that should be appointed, that’s the most used version of _D“ is when 
a = 0, so 

D ° m = Taho) (!) I {t ~ xr ~°" f (x)Jx ■ 

In 1967, M. Caputo presented a new definition of a fractional derivative called Caputo’s Fractional 
Derivative which is a modification of the Riemann-Liouville Fractional Derivative [6]. However, the 
definition of the Caputo’s Fractional Derivative uses almost the same notations of the definition of 
the Riemann-Liouville Fractional Derivative. 


Definition 2.3. Let a G M + and n G N such that n — 1 < a < n, then the Caputo’s Fractional 
Derivative of order a is defined by 


D a a f (t) = 


r(n-a ) J a (;t - X) 


f {n) jx) 

a— n+1 


.dx. 


( 2 , 6 ) 


In addition, if a = 0 in Eq. (2.6), then the most used version of the Caputo’s Fractional Derivative 
is obtained, i.e. 


D a f (t) = 


f(n) 


X 


r(n-a) J 0 ( ; t - x) 


a—n+1 


.dx. 


(2.7) 


3. Approximating the Fractional Derivative of order 1 < ot < 2 using FFDM 


For the fractional derivative, the Caputo’s definition is chosen. Now, according to Eq. (2.7), if 
1 < a < 2, then the Caputo’s Fractional Derivative is defined by 


D a y (f) = 


y" (x) 


r (2- a ) J 0 (t - x)‘ 


— ■dx 


(3.1) 


for t > 0 and a G M + . 

By applying integral by parts on the right hand side of Eq. (3.1), we get 

1 


y' (0) t 2 a + / (t — x) 2 a y "' (x) .dx ) . 


D a y (t) 


(2 - a) r (2 - a ) 


(3.2) 
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Last integral in Eq. (3.2) can be approximated using the Composite Trapezoidal Rule, Eq. (2.5), as 
follows 


(t ~ xf a y" ( x) Ax » ^ t 2 a y" (0) + 2 ^ (t - Xj ) 2 a y" ( Xj ) + (t - bf a y" (6)], (3.3) 

3 = 1 


n—1 


with h — {b — a)/n and Xj = a + jh for each j — 0,1,..., n — 1. 
Now, by substituting Eq. (3.3) in Eq. (3.2), we obtain, 


D a y(t ) 


(2 - a) r (2 - a ) 


y" (0) t 


2—a 


h ; 
2 ' 


n—1 


t 2 a y" (o) + 2 - x jf a y" ( x j) + ( t ~ h f a y" ( & )] 

3 = 1 


(3.4) 


Here, y "and y "can be approximated at specihc points. Forward finite difference formula for the 2 nd 

i to ap 

V (0) 


order derivative, Eq. (2.1), is used to approximate y (0) such that 

y (0) - 2 y (, h ) + y (2 h) 


h 2 


(3.5) 


for small values of h. 

To approximate y (0), also forward finite difference formula for the 3 rd derivative, Eq. (2.3), is used 
as follows 

tr\\ _ -y (°) + 3 y( h ) - 3 y (2 h) + y(3h) 
y (0) ~- (3.6) 


h 3 


for small values of h. 


While the central hnite difference formula for the 3 rd derivative, Eq. (2.2), is used to approximate 
y {xj) in the sum term of Eq. (|3.4) as 


y" {xj) 


-y [xj — 2 h) + 2 y( X j — h) — 2y {xj + h) + y(xj + 2 h) 

2 h 3 


(3.7) 


for small values of h and for j = 0,1,..., n — 1. 

Finally, y ( b ) can be approximated using backward hnite difference formula for the 3 rd derivative, 


Eq. (2.4), as follows 


y" (6) 


y ( b ) — 3 y(b— h ) + 3 y [b — 2 h) — y(b — 3 h) 

h 3 


(3.8) 


for small values of h. 

By substituting Eq.s (3.5), (3.6), (3.7) and (3.8) in Eq. (3.4), we get 


D a y W « ^-r(A + \(B + 2 C + D)), 

(2 — a) T (2 — a ) 2 


(3.9) 


, _ y(0) - 2y(h) + y(2h) t2 _ a 
A ~ ¥ * ’ 

-y (0) + 3 y(h) - 3 y (2 h) + y(3h) 2 _ a 
B = - ¥ -‘ ’ 


where 
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and 


n— 1 


j =i 



2 h) + 2 y(xj 


h ) — 2 y (xj + h ) + + 2/r) 

~2/D 


(t - x,) 2 -“ 


5 


D = 


y (b) -3 y(b-h) + 3 y (b - 2h) - y(b - 3h) 

h 3 


(t-6) 


2—a 


We observe that, according to Eq. (3.9), the approximation of D a (y(t )) depends on the value 


of h such that if the value of h is decreased, then the approximate result of D a y (t) will be more 


accurate. Moreover, Eq. (3.9) is the purposed equation that needs an algorithm to approximate the 
fractional derivative of order 1 < a < 2 of a given function y(t) which is clear. Also, the value of h 


according to the Composite Trapezoidal Rule, Eq. (2.5), is (b — a)/n So; 


n = 


h 


Now, we can find D a y (t) approximately by taking increasingly n’s, that’s we may take n={10, 100, 
1000, 10000}. Also, we should note that we use a = 0 and 6 = 1. 


4. Numerical Solutions for Linear Fractional Differential Equations of order 1 < a < 2 
using FFDM 

To show the efficiency of the FFDM, we will approximate the solutions for some linear fractional 
differential equations of order 1 < a < 2 using a prepared program in Mathematica (Ver. 9.0) 
and compare them with the exact solutions and with other approximate solutions obtained by other 
solvers and methods. As we mentioned before, Eq. (3.9) is considered as an approximation for 


D a y (f); so whenever we find D a y (t) in any linear fractional differential equation of order 1 < a < 2 


of the form of Eq. (1.1), we may replace it by that approximation. 

Example 4.1 (|22j). Consider the following linear fractional differential equation: 


D 1 - 5 y(t)=t 1 - 5 .y(t)+4 ] f^-t 3 -\ 
with the following initial conditions 

2 /( 0 ) = 0 . 

Note that, the exact solution for this problem is y ( t ) = t 2 . 

By taking h = 0.1, Table [l] illustrates the concluded results using FFDM by computing the 
approximate solutions for this problem at specific points of t. 

Example 4.2 (pf, 27]). Consider the following linear fractional differential equation: 

D 2 y (t) + D ' y (t) + y (t) = t 2 + 2 + 4^-j, 
with the following initial condition 

2 /( 0 ) = 0 . 

Note that, the exact solution for this problem is y ( t ) = t 2 . 

Here, we take h = 0.1, Table [2] illustrates the concluded results using FFDM by computing the 
approximate solutions for this problem at specific points of t. 
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Table 1: Numerical solutions for Example 4.1 using FFDM, when h = 0.1 


T 

App. 

Exact 

Absolute Error 

0.1 

0.00999999999999998 

0.01 

2.08167E-17 

0.2 

0.03999999999999981 

0.04 

2.01228E-16 

0.3 

0.08999999999999946 

0.09 

5.96745E-16 

0.4 

0.15999999999999906 

0.16 

9.99201E-16 

0.5 

0.24999999999999858 

0.25 

1.9984E-15 

0.6 

0.35999999999999790 

0.36 

2.9976E-15 

0.7 

0.48999999999999705 

0.49 

2.9976E-15 

0.8 

0.63999999999999600 

0.64 

3.9968E-15 

0.9 

0.80999999999999470 

0.81 

6.10623E-15 

1.0 

0.99999999999999340 

1.00 

6.99441E-15 


Table 2: Numerical solutions for Example 4.2 using FFDM, when h = 0.1 


T 

App. 

Exact 

Absolute Error 

0.1 

0.01000000000000325 

0.01 

3.20056E-15 

0.2 

0.04000000000000937 

0.04 

9.29812E-15 

0.3 

0.09000000000001489 

0.09 

1.48076E-14 

0.4 

0.16000000000001902 

0.16 

1.89848E-14 

0.5 

0.25000000000002354 

0.25 

2.29816E-14 

0.6 

0.36000000000002824 

0.36 

2.80331E-14 

0.7 

0.49000000000003320 

0.49 

3.30291E-14 

0.8 

0.64000000000003820 

0.64 

3.79696E-14 

0.9 

0.81000000000004320 

0.81 

4.29656E-14 

1.0 

1.00000000000004860 

1.00 

3.9968E-14 


Example 4.3 ([13]). Consider the following linear fractional differential equation: 

D 2 y (t) + D L5 y (t) + y(t) = l + t, 
with the following initial condition 

1 /( 0 ) = 1 - 

Regarding to [13] . the author’s develop a new method called Generalized Taylor Collocation Method 
(GTCM) to give approximate solutions for linear fractional differential equations with variable co¬ 
efficients. They found that this method is easy to write in computer codes and can be used as an 
effective method for obtaining analytic and approximate solutions for fractional differential equations. 
Moreover, the previous problem was solved within their paper by their method. They introduce the 
fourth order approximation solution for that problem which is as 

y(t) = i + t. 

However, our FFDM is also applied to solve this problem and compared with GTCM according to 
their approximate solutions, by taking h = 0.1, Figure [l] shows that the concluded results using 
FFDM and the results using GTCM coincide. 
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Figure 1: The comparison between the FFDM solutions (where h = 0.1) and the approximate solutions which are 
obtained by GTCM, according to Example |4.3 


5. Conclusions 

The major goal of this paper is to find accurate approximate solutions for linear fractional dif¬ 
ferential equations of order 1 < a < 2. Hence, we carry out this goal by preparing a new method 
called Fractional Finite Difference Method (FFDM). In this paper, we discussed and illustrated the 
numerical solutions of linear fractional differential equations using FFDM. Also, the efficiency of 
FFDM was illustrated by solving some examples of linear fractional differential equations of order 
1 < a < 2 . 

Successfully, FFDM was applied to approximate the linear fractional differential equations of 
order 1 < a < 2. All ideas were illustrated to be efficient in applying the proposed technique 
to several linear fractional differential equations of order 1 < a < 2. We found that our method 
is powerful and efficient in finding numerical solutions for those equations. Moreover, the FFDM 
solutions demonstrate excellent approximations in comparison with the exact solutions and with 
other methods and solvers through the applicable domain. 
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